EP playground – 06 categ dummy

library(tidyverse)
library(patchwork)
library(emmeans)
library(simglm)
library(ggforce)

Simulate data

set.seed(3119) 

sim_arguments <- list(
  formula = y ~ 1 + hours + motivation + study + method,
  fixed = list(hours = list(var_type = 'ordinal', levels = 0:15),
               motivation = list(var_type = 'continuous', mean = 0, sd = 1),
               study = list(var_type = 'factor', 
                            levels = c('alone', 'others'),
                            prob = c(0.53, 0.47)),
               method = list(var_type = 'factor', 
                            levels = c('read', 'summarise', 'self-test'),
                            prob = c(0.3, 0.4, 0.3))),
  error = list(variance = 20),
  sample_size = 250,
  reg_weights = c(0.6, 1.4, 1.5, 6, 6, 2)
)

df3 <- simulate_fixed(data = NULL, sim_arguments) %>%
  simulate_error(sim_arguments) %>%
  generate_response(sim_arguments)

test_study3 <- df3 %>%
  dplyr::select(y, hours, motivation, study, method) %>%
  mutate(
    ID = paste("ID", 101:350, sep = ""),
    score = round(y+abs(min(y))),
    motivation = round(motivation, 2),
    study = factor(study),
    method = factor(method)
  ) %>%
  dplyr::select(ID, score, hours, motivation, study, method)

study (two levels)

Group means

test_study3 |>
  group_by(study) |>
  summarise(
    avg_score = round(mean(score), 2)
  )

mean_alone <- filter(test_study3, study == 'alone')$score |> mean()
mean_others <- filter(test_study3, study == 'others')$score |> mean()

Plot study data

test_study3 |>
  ggplot(aes(x = study, y = score, fill = study, colour = study)) +
  geom_violin(alpha = 0.5) +
  geom_sina(alpha = 0.5) +
  theme(legend.position = 'none') +
  stat_summary(colour = 'black', fun = mean, geom = 'point', size = 2) +
  geom_segment(aes(x = 'alone', xend = 'others', y = mean_alone, yend = mean_others), colour = 'black') +
  NULL

Plot on xy plane, dummy-coded.

contrasts(test_study3$study)

       others
alone       0
others      1

alone = ref level.

xlim_lower <- -2.2
xlim_upper <-  2.2
ylim_lower <- -25
ylim_upper <-  55

test_study3 |>
  mutate(study_num = ifelse(study == 'alone', 0, 1)) |>
  ggplot(aes(x = study_num, y = score, fill = study, colour = study)) +
  geom_jitter(alpha = 0.5, width = 0.1) +
  scale_x_continuous(limits = c(xlim_lower, xlim_upper), expand = c(0, 0)) +
  scale_y_continuous(limits = c(ylim_lower, ylim_upper), expand = c(0, 0)) +
  geom_segment(aes(x = xlim_lower, xend = xlim_upper, y = 0, yend = 0), arrow = arrow(ends = 'both', length = unit(12, 'pt')), colour = 'black') +
  geom_segment(aes(x = 0, xend = 0, y = ylim_lower, yend = ylim_upper), arrow = arrow(ends = 'both', length = unit(12, 'pt')), colour = 'black') +
  geom_segment(aes(x = 0, xend = 1, y = mean_alone, yend = mean_others), colour = 'black') +
  geom_segment(aes(x = 0, xend = 1, y = mean_alone, yend = mean_alone), colour = 'red', linewidth = 2) +
  geom_segment(aes(x = 1, xend = 1, y = mean_alone, yend = mean_others), colour = 'red', linewidth = 2) +
  NULL

Add note that the manual stuff is so you understand how coding works. But you’ll never have to do it this laboriously.

Exact same outcome:

• manual ifelse
• contr.treatment()

Model score ~ study

mod1 <- lm(score ~ study, data = test_study3)
summary(mod1)

Call:
lm(formula = score ~ study, data = test_study3)

Residuals:
     Min       1Q   Median       3Q      Max 
-22.7902  -4.7902   0.2098   5.4766  18.2098 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  22.7902     0.6603   34.52  < 2e-16 ***
studyothers   4.7332     1.0093    4.69 4.53e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.896 on 248 degrees of freedom
Multiple R-squared:  0.08145,   Adjusted R-squared:  0.07775 
F-statistic: 21.99 on 1 and 248 DF,  p-value: 4.525e-06

t-test

t.test(score ~ study, var.equal = T, data = test_study3)

    Two Sample t-test

data:  score by study
t = -4.6895, df = 248, p-value = 4.525e-06
alternative hypothesis: true difference in means between group alone and group others is not equal to 0
95 percent confidence interval:
 -6.721058 -2.745251
sample estimates:
 mean in group alone mean in group others 
            22.79021             27.52336 

emmeans

(mod1_emm <- emmeans(mod1, ~study))

 study  emmean    SE  df lower.CL upper.CL
 alone    22.8 0.660 248     21.5     24.1
 others   27.5 0.763 248     26.0     29.0

Confidence level used: 0.95 

(p_mod1_emm <- mod1_emm |> plot() +
  coord_flip())

overlay emmeans on data plot

mod1_emm_tib <- mod1_emm |>
  as_tibble() |>
  rename(score = emmean)

test_study3 |>
  ggplot(aes(x = study, y = score, fill = study, colour = study)) +
  geom_violin(alpha = 0.5) +
  geom_jitter(alpha = 0.5, width = 0.2) +
  theme(legend.position = 'none') +
  geom_errorbar(data = mod1_emm_tib, aes(ymin = lower.CL, ymax = upper.CL), colour = 'black', width = 0.2) +
  geom_point(data = mod1_emm_tib, colour = 'black', size = 3) +
  NULL

pairwise comparisons

levels(test_study3$study)

[1] "alone"  "others"

For each comparison we want to make, assign 1 to the level we’re interested in, and assign -1 to the level you want to compare it to. If there are any other levels in there, assign them 0.

mod1_comparisons <- list(
  'alone vs. others' = c(
    1,   # alone
    -1   # others
  )
)

contrast(
  object = mod1_emm,
  method = mod1_comparisons
)

 contrast         estimate   SE  df t.ratio p.value
 alone vs. others    -4.73 1.01 248  -4.690  <.0001

Basically reconstructs the model’s estimates. Not so interesting.

Change ref level: coefs are diff but emmeans are same

Original contrasts:

contrasts(test_study3$study)

       others
alone       0
others      1

Now, let’s flip it:

contrasts(test_study3$study) <- contr.treatment(
  levels(test_study3$study), 
  base = 2)
contrasts(test_study3$study)

       alone
alone      1
others     0

plot flipped data

xlim_lower <- -2.2
xlim_upper <-  2.2
ylim_lower <- -25
ylim_upper <-  55

test_study3 |>
  mutate(study_num = ifelse(study == 'alone', 1, 0)) |>
  ggplot(aes(x = study_num, y = score, fill = study, colour = study)) +
  geom_jitter(alpha = 0.5, width = 0.1) +
  scale_x_continuous(limits = c(xlim_lower, xlim_upper), expand = c(0, 0)) +
  scale_y_continuous(limits = c(ylim_lower, ylim_upper), expand = c(0, 0)) +
  geom_segment(aes(x = xlim_lower, xend = xlim_upper, y = 0, yend = 0), arrow = arrow(ends = 'both', length = unit(12, 'pt')), colour = 'black') +
  geom_segment(aes(x = 0, xend = 0, y = ylim_lower, yend = ylim_upper), arrow = arrow(ends = 'both', length = unit(12, 'pt')), colour = 'black') +
  geom_segment(aes(x = 0, xend = 1, y = mean_others, yend = mean_alone), colour = 'black') +
  geom_segment(aes(x = 0, xend = 1, y = mean_others, yend = mean_others), colour = 'red', linewidth = 2) +
  geom_segment(aes(x = 1, xend = 1, y = mean_others, yend = mean_alone), colour = 'red', linewidth = 2) +
  NULL

fit model (different)

mod1b <- lm(score ~ study, data = test_study3)
coef(mod1b)

(Intercept)  studyalone 
  27.523364   -4.733155 

• What is intercept? (mean when study = others)
• What is studyalone? (diff between levels: mean score when study = alone minus mean score when study = others)

Compare this to the original coefs from mod1:

coef(mod1)

(Intercept) studyothers 
  22.790210    4.733155 

• What is intercept? (mean when study = alone)
• What is studyothers? (diff between levels: mean score when study = others minus mean score when study = alone)

get emmeans (same)

(mod1b_emm <- emmeans(mod1b, ~study))

 study  emmean    SE  df lower.CL upper.CL
 alone    22.8 0.660 248     21.5     24.1
 others   27.5 0.763 248     26.0     29.0

Confidence level used: 0.95 

p_mod1b_emm <- mod1b_emm |> plot() +
  coord_flip() +
  ggtitle('With reference level = others')

p_mod1_emm <- p_mod1_emm + ggtitle('With reference level = alone')

p_mod1b_emm + p_mod1_emm

Identical!

method (three levels)

Group means

test_study3 |>
  group_by(method) |>
  summarise(
    avg_score = round(mean(score), 2)
  )

mean_read <- filter(test_study3, method == 'read')$score |> mean()
mean_self <- filter(test_study3, method == 'self-test')$score |> mean()
mean_summ <- filter(test_study3, method == 'summarise')$score |> mean()

Plot data

test_study3 |>
  ggplot(aes(x = method, y = score, fill = method, colour = method)) +
  geom_violin(alpha = 0.5) +
  geom_jitter(alpha = 0.5, width = 0.2) +
  theme(legend.position = 'none') +
  NULL

Goal: Compare the means of each group to one another. That was easy when we had only two groups: one straight line can go from mean to mean. But now one straight line won’t cut it anymore.

The smallest number of lines that will connect all the means = 2. Instead of estimating one line, we will estimate two.

And in general, if we have \(n\) groups, we will be estimating \(n-1\) lines.

The lines that we’re going to be estimating now:

test_study3 |>
  ggplot(aes(x = method, y = score, fill = method, colour = method)) +
  geom_violin(alpha = 0.5) +
  geom_jitter(alpha = 0.5, width = 0.2) +
  theme(legend.position = 'none') +
  stat_summary(colour = 'black', fun = mean, geom = 'point', size = 2) +
  # line from read to self-test:
  geom_segment(colour = 'black', aes(x = 'read', xend = 'self-test', y = mean_read, yend = mean_self)) +
  # line from read to summarise:
  geom_segment(colour = 'black', aes(x = 'read', xend = 'summarise', y = mean_read, yend = mean_summ)) +
  NULL

But I’m lying a little bit. This looks like going to go from 0 to 1 in one line, and 0 to 2 in the other line. But that’s not true – this is just the intuition about what groups we are comparing. Really what we’re going to do is this:

Plot dummy-coded method in xy space

Two panels: the first will be the difference between read and self-test (the first variable), the second will be difference between read and summarise (the second variable).

xlim_lower <- -2.2
xlim_upper <-  2.2
ylim_lower <- -25
ylim_upper <-  55

p1 <- test_study3 |> 
  filter(method %in% c('read', 'self-test')) |>
  mutate(method_num = ifelse(method == 'read', 0, 1)) |>
  ggplot(aes(x = method_num, y = score, fill = method, colour = method)) +
  geom_jitter(alpha = 0.5, width = 0.1) +
  scale_x_continuous(limits = c(xlim_lower, xlim_upper), expand = c(0, 0)) +
  scale_y_continuous(limits = c(ylim_lower, ylim_upper), expand = c(0, 0)) +
  geom_segment(aes(x = xlim_lower, xend = xlim_upper, y = 0, yend = 0), arrow = arrow(ends = 'both', length = unit(12, 'pt')), colour = 'black') +
  geom_segment(aes(x = 0, xend = 0, y = ylim_lower, yend = ylim_upper), arrow = arrow(ends = 'both', length = unit(12, 'pt')), colour = 'black') +
  geom_segment(aes(x = 0, xend = 1, y = mean_read, yend = mean_self), colour = 'black') +
  geom_segment(aes(x = 0, xend = 1, y = mean_read, yend = mean_read), colour = 'red', linewidth = 2) +
  geom_segment(aes(x = 1, xend = 1, y = mean_read, yend = mean_self), colour = 'red', linewidth = 2) +
  labs(
    title = 'First line: read vs. self-test'
  ) +
  NULL

p2 <- test_study3 |> 
  filter(method %in% c('read', 'summarise')) |>
  mutate(method_num = ifelse(method == 'read', 0, 1)) |>
  ggplot(aes(x = method_num, y = score, fill = method, colour = method)) +
  geom_jitter(alpha = 0.5, width = 0.1) +
  scale_x_continuous(limits = c(xlim_lower, xlim_upper), expand = c(0, 0)) +
  scale_y_continuous(limits = c(ylim_lower, ylim_upper), expand = c(0, 0)) +
  geom_segment(aes(x = xlim_lower, xend = xlim_upper, y = 0, yend = 0), arrow = arrow(ends = 'both', length = unit(12, 'pt')), colour = 'black') +
  geom_segment(aes(x = 0, xend = 0, y = ylim_lower, yend = ylim_upper), arrow = arrow(ends = 'both', length = unit(12, 'pt')), colour = 'black') +
  geom_segment(aes(x = 0, xend = 1, y = mean_read, yend = mean_summ), colour = 'black') +
  geom_segment(aes(x = 0, xend = 1, y = mean_read, yend = mean_read), colour = 'red', linewidth = 2) +
  geom_segment(aes(x = 1, xend = 1, y = mean_read, yend = mean_summ), colour = 'red', linewidth = 2) +
  labs(
    title = 'Second line: read vs. summarise'
  ) +  
  NULL

p1 + p2

Dummy-code method variable

Manual dummy-coding:

dummyCoded <- test_study3 %>%
  select(ID, score, method) %>%
  mutate(
    method1 = ifelse(method == "self-test", 1, 0),
    method2 = ifelse(method == "summarise", 1, 0))

dummyCoded

contrasts(test_study3$method)

          self-test summarise
read              0         0
self-test         1         0
summarise         0         1

Model score ~ method

mod2 <- lm(score ~ method, data = dummyCoded)
summary(mod2)

Call:
lm(formula = score ~ method, data = dummyCoded)

Residuals:
     Min       1Q   Median       3Q      Max 
-23.4138  -5.3593  -0.1959   5.7496  17.8041 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)      23.4138     0.8662  27.031   <2e-16 ***
methodself-test   4.1620     1.3188   3.156   0.0018 ** 
methodsummarise   0.7821     1.1930   0.656   0.5127    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.079 on 247 degrees of freedom
Multiple R-squared:  0.04224,   Adjusted R-squared:  0.03448 
F-statistic: 5.447 on 2 and 247 DF,  p-value: 0.004845

t-tests

(Nearly) equivalent t-tests (different bc diff degrees of freedom, bc filtering dataset. but close).

dC1 <- test_study3 |> 
  filter(method %in% c('read', 'self-test')) |>
  mutate(method = factor(method, levels = c('read', 'self-test')))
contrasts(dC1$method) <- contr.treatment(2)
contrasts(dC1$method)

          2
read      0
self-test 1

t.test(score ~ method, var.equal = T, data = dC1)

    Two Sample t-test

data:  score by method
t = -3.1354, df = 151, p-value = 0.002063
alternative hypothesis: true difference in means between group read and group self-test is not equal to 0
95 percent confidence interval:
 -6.784657 -1.539272
sample estimates:
     mean in group read mean in group self-test 
               23.41379                27.57576 

dC2 <- test_study3 |> 
  filter(method %in% c('read', 'summarise')) |>
  mutate(method = factor(method, levels = c('read', 'summarise')))
contrasts(dC2$method) <- contr.treatment(2)
contrasts(dC2$method)

          2
read      0
summarise 1

t.test(score ~ method, var.equal = T, data = dC2)

    Two Sample t-test

data:  score by method
t = -0.66342, df = 182, p-value = 0.5079
alternative hypothesis: true difference in means between group read and group summarise is not equal to 0
95 percent confidence interval:
 -3.108082  1.543915
sample estimates:
     mean in group read mean in group summarise 
               23.41379                24.19588 

emmeans

Before, I said that the smallest number of lines we need to compare three groups is two. But this actually means that one of the possible comparisons is left out.

(??? above the missing line between self-test and summarise)

Options:

• make something else ref level. Then we get that comparison. But then we’re missing a different third comparison.
• get emmeans, and once we’re in the outcome space, we can compute differences between whatever we want! yay!

(mod2_emm <- emmeans(mod2, ~method))

 method    emmean    SE  df lower.CL upper.CL
 read        23.4 0.866 247     21.7     25.1
 self-test   27.6 0.994 247     25.6     29.5
 summarise   24.2 0.820 247     22.6     25.8

Confidence level used: 0.95 

mod2_emm |> plot() +
  coord_flip()

overlay emmeans on data plot

mod2_emm_tib <- mod2_emm |>
  as_tibble() |>
  rename(score = emmean)

test_study3 |>
  ggplot(aes(x = method, y = score, fill = method, colour = method)) +
  geom_violin(alpha = 0.5) +
  geom_jitter(alpha = 0.5, width = 0.2) +
  theme(legend.position = 'none') +
  geom_errorbar(data = mod2_emm_tib, aes(ymin = lower.CL, ymax = upper.CL), colour = 'black', width = 0.2) +
  geom_point(data = mod2_emm_tib, colour = 'black', size = 3) +
  NULL

pairwise comparisons

levels(test_study3$method)

[1] "read"      "self-test" "summarise"

Assign 1 to the level we are interested in, and -1 to the baseline we want to compare it to.

mod2_comparisons <- list(
  'Self-test vs. Read' = c(-1, 1, 0),  # read self-test summarise, in that order
  'Summarise vs. Read' = c(-1, 0, 1),
  'Self-test vs. Summarise' = c(0, 1, -1),
  'mean(Self-test, summarise) vs. Read' = c(-1, 0.5, 0.5)  # weights must sum to 0 -- pools self-test and summarise tog and pits that mean against mean of Read
)

1 - 1

mean(Self-test) - mean(Read) = positive, bc self-test is larger than read

(mod2_comparisons_test <- contrast(mod2_emm, mod2_comparisons))

 contrast                            estimate   SE  df t.ratio p.value
 Self-test vs. Read                     4.162 1.32 247   3.156  0.0018
 Summarise vs. Read                     0.782 1.19 247   0.656  0.5127
 Self-test vs. Summarise                3.380 1.29 247   2.622  0.0093
 mean(Self-test, summarise) vs. Read    2.472 1.08 247   2.290  0.0229

contrast() gives us the p-values for each difference. We can also get the 95% CIs of each difference by giving the outcome of contrast() into confint().

confint(mod2_comparisons_test)

 contrast                            estimate   SE  df lower.CL upper.CL
 Self-test vs. Read                     4.162 1.32 247    1.564     6.76
 Summarise vs. Read                     0.782 1.19 247   -1.568     3.13
 Self-test vs. Summarise                3.380 1.29 247    0.841     5.92
 mean(Self-test, summarise) vs. Read    2.472 1.08 247    0.345     4.60

Confidence level used: 0.95